1 Objective

The goal of this session is to get a clean macula densa dataset to work with.

2 Load Packages

options(future.globals.maxSize = 74 * 1024^3) # 55 GB
getOption("future.globals.maxSize") #59055800320
## [1] 79456894976

3 Load Dataset

SO1 <- LoadSeuratRds(here("outputs", "so_merge_raw.rds"))


head(SO1@meta.data)

4 Subset out non MD cells.

SO2 <- subset(SO1, idents = c("4", "5", "6", "7"), invert = TRUE)

SO3 <- subset(SO2, features = grep("^mt-|^Gm|Rik|^Rp", rownames(SO2), invert = TRUE, value = TRUE))

DimPlot(SO2)

SO3 <- SCTransform(SO3) %>%
    RunPCA() %>%
    FindNeighbors(dims = 1:15) %>%
    FindClusters(resolution = 2) %>%
    RunUMAP(dims = 1:15)
## Running SCTransform on assay: RNA
## Warning: The `slot` argument of `GetAssayData()` is deprecated as of SeuratObject 5.0.0.
## ℹ Please use the `layer` argument instead.
## ℹ The deprecated feature was likely used in the Seurat package.
##   Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.
## Calculating cell attributes from input UMI matrix: log_umi
## Variance stabilizing transformation of count matrix of size 14658 by 11529
## Model formula is y ~ log_umi
## Get Negative Binomial regression parameters per gene
## Using 2000 genes, 5000 cells
## There are 2 estimated thetas smaller than 1e-07 - will be set to 1e-07
## Found 302 outliers - those will be ignored in fitting/regularization step
## Second step: Get residuals using fitted parameters for 14658 genes
## Computing corrected count matrix for 14658 genes
## Calculating gene attributes
## Wall clock passed: Time difference of 15.0815 secs
## Determine variable features
## Centering data matrix
## Place corrected count matrix in counts slot
## Warning: The `slot` argument of `SetAssayData()` is deprecated as of SeuratObject 5.0.0.
## ℹ Please use the `layer` argument instead.
## ℹ The deprecated feature was likely used in the Seurat package.
##   Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: Different cells and/or features from existing assay SCT
## Set default assay to SCT
## PC_ 1 
## Positive:  Umod, Egf, Tmem52b, Sult1d1, Sostdc1, Fabp3, Prdx5, Foxq1, Klk1, Cox6c 
##     Wfdc15b, Krt7, Ly6a, Ckb, Slc25a5, Cox5b, Wfdc2, Atp5g1, Cldn19, Ggt1 
##     Ndufa4, Cox4i1, Atp1a1, Chchd10, Atp5b, Cox8a, Gadd45g, Cox6b1, Cox7a2, Uqcr11 
## Negative:  Pappa2, Zfand5, Aard, Robo2, Pde10a, Wwc2, Itga4, Nadk2, S100g, Nos1 
##     Ramp3, Neat1, Sgms2, Ptgs2, Col4a4, Irx1, Col4a3, Tmem158, Mir6236, Ranbp3l 
##     Itprid2, Bmp3, Cdkn1c, Camk2d, Dctd, Mcub, Sdc4, Etnk1, Peg3, Zbtb20 
## PC_ 2 
## Positive:  Fos, Junb, Jun, Egr1, Btg2, Hspa1b, Hspa1a, Atf3, Zfp36, Ier2 
##     Fosb, Socs3, Klf6, Klf2, Jund, Dusp1, Dnajb1, Rhob, Gadd45g, Tsc22d1 
##     Gadd45b, Klf4, H3f3b, Ubc, Cebpd, Actb, H1f2, H2bc4, Ppp1r15a, Hspb1 
## Negative:  Pappa2, Egf, Slc12a1, CT010467.1, Sfrp1, Mir6236, Etnk1, Umod, Wnk1, Sec14l1 
##     Robo2, Hsp90b1, Pde10a, Car15, Atp1b1, Col4a4, Aard, Wwc2, Mal, Itprid2 
##     Mgst1, Aebp1, Col4a3, Bmp3, Trim2, Nos1, Sgms2, Zbtb20, Pth1r, Tfcp2l1 
## PC_ 3 
## Positive:  Fth1, Ubb, Ftl1, Ldhb, Fxyd2, Car15, Prdx1, Mt1, Eif1, Cd63 
##     Cd9, Mpc2, Mgst1, Clu, Bsg, Tmem213, Aldoa, Mdh1, Ramp3, S100a1 
##     Selenow, Spp1, Tmem176b, Tspo, H3f3b, Wfdc2, Atpif1, Lgals3, Tmem59, Tmbim6 
## Negative:  Egf, Mir6236, CT010467.1, Umod, Tmem52b, Neat1, Nme7, Malat1, Slc12a1, Kcnq1ot1 
##     Lars2, Etnk1, Wnk1, Dst, Sult1d1, Foxq1, Slc5a3, Atrx, Atp1b1, Pnisr 
##     Fos, Syne2, Gls, Ivns1abp, Hnrnpa2b1, Srrm2, Chka, Col4a4, Paxbp1, Zbtb20 
## PC_ 4 
## Positive:  Mt1, Apoe, Mt2, Aebp1, Sostdc1, Gpx6, Fxyd2, Neat1, Ptger3, Ckb 
##     Fgf9, Egfl6, Ivns1abp, Defb1, Mfsd4a, Car4, Fkbp11, Chchd10, Plet1, Hspa1b 
##     Atp5md, Igfbp5, Cox6c, Tmem213, Atp5k, CT010467.1, Chka, Mgst3, Lpl, Atp1a1 
## Negative:  Pappa2, Aard, Tmem52b, Umod, Sult1d1, Tmem158, Egf, Foxq1, Mcub, Ptgs2 
##     Dctd, Wwc2, Iyd, Ramp3, S100g, Ptprz1, Cd9, Car15, Hsp90b1, Tmsb4x 
##     Defb42, Cdkn1c, Wnk1, Pth1r, Tagln2, Lcn2, Il1f6, Clu, Ctsc, Bmp2 
## PC_ 5 
## Positive:  Hspa1b, Hspa1a, Umod, Egf, Jun, Cldn10, Pappa2, Car15, Klk1, Fos 
##     Sfrp1, Fth1, Ier3, Clcnkb, Atp1a1, Slc12a1, Cited2, Hspa8, Pik3r1, Cdkn1c 
##     Robo2, Kng2, Aard, Gsta4, Ptger3, Fau, Irx1, Ftl1, Uba52-ps, Tmem37 
## Negative:  S100g, Actb, Tmem52b, Igfbp7, Cryab, S100a6, Tacstd2, Tmsb10, Tmsb4x, Atf3 
##     Lcn2, Sult1d1, Foxq1, Crip1, Egr1, Ppia, Atp5md, Lgals1, Mir6236, Gadd45b 
##     Klf4, Ifitm3, Ccnd1, Ccn1, Hbegf, Serf2, Spp1, Nupr1, Cdkn1a, Krt7
## Computing nearest neighbor graph
## Computing SNN
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 11529
## Number of edges: 359040
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7101
## Number of communities: 25
## Elapsed time: 1 seconds
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
## 13:27:36 UMAP embedding parameters a = 0.9922 b = 1.112
## 13:27:36 Read 11529 rows and found 15 numeric columns
## 13:27:36 Using Annoy for neighbor search, n_neighbors = 30
## 13:27:36 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 13:27:37 Writing NN index file to temp file /var/folders/7n/x74qctp91rng390gx0z9hmd80000gn/T//RtmpIFjmQt/file54ff5c8bcfcc
## 13:27:37 Searching Annoy index using 1 thread, search_k = 3000
## 13:27:39 Annoy recall = 100%
## 13:27:39 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 13:27:40 Initializing from normalized Laplacian + noise (using RSpectra)
## 13:27:40 Commencing optimization for 200 epochs, with 484894 positive edges
## 13:27:40 Using rng type: pcg
## 13:27:42 Optimization finished
DimPlot(SO3)

DimPlot(SO3,split.by = "treatment")

markers.to.plot1 <- c("Lrp2",         # PT
                      "Slc5a12",      # PT-S1
                      "Slc13a3",      # PT-S2
                      "Slc16a9",      # PT-S3
                      "Havcr1",       # Injured PT
                      "Epha7",        # dTL
                      "Cryab",        # dTL
                      "Cdh13",        # dTL1
                      "Slc14a2",      # dTL2
                      "Slc12a1",      # TAL
                      "Umod",         # TAL, DCT1
                      "Egf",          # TAL, DCT1,
                      "Cldn10",       # TAL
                      "Cldn16",       # TAL
                      "Nos1",         # MD
                      "Slc12a3",      # DCT
                      "Pvalb",        # DCT1
                      "Slc8a1",       # DCT2, CNT
                      "Aqp2",         # PC
                      "Slc4a1",       # IC-A
                      "Slc26a4",      # IC-B
                      "Nphs1",        # Podo
                      "Ncam1",        # PEC
                      "Flt1",         # Endo
                      "Emcn",         # Glom Endo
                      "Kdr",          # Capillary Endo
                      "Pdgfrb",       # Perivascular
                      "Pdgfra",       # Fib
                      "Piezo2",       # Mesangial
                      "Acta2",        # Mural
                      "Ptprc",        # Immune
                      "Cd74",         # Macrophage
                      "Skap1",        # B/T Cells 
                      "Upk1b",        # Uro
                      "Top2a"         # Proliferation
)

DotPlot(SO3,
features = markers.to.plot1,
dot.scale = 8,
dot.min = 0,
scale = FALSE,
scale.max = 100,
scale.min = 0,
col.min = -2.5,
col.max = 2.5)+
coord_flip()
## Warning: Could not find Slc5a12 in the default search locations, found in 'RNA'
## assay instead
## Warning: Could not find Kdr in the default search locations, found in 'RNA'
## assay instead
## Warning: The following requested variables were not found: Slc13a3, Havcr1,
## Pdgfrb, Pdgfra, Piezo2, Upk1b

head(SO3)
SO4 <- subset(SO3, idents = c("23"), invert = TRUE)

DimPlot(SO4)

SO4 <- SCTransform(SO4) %>%
    RunPCA() %>%
    FindNeighbors(dims = 1:16) %>%
    FindClusters(resolution = 2) %>%
    RunUMAP(dims = 1:16)
## Running SCTransform on assay: RNA
## vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.
## Calculating cell attributes from input UMI matrix: log_umi
## Variance stabilizing transformation of count matrix of size 14577 by 11448
## Model formula is y ~ log_umi
## Get Negative Binomial regression parameters per gene
## Using 2000 genes, 5000 cells
## There are 2 estimated thetas smaller than 1e-07 - will be set to 1e-07
## Found 292 outliers - those will be ignored in fitting/regularization step
## Second step: Get residuals using fitted parameters for 14577 genes
## Computing corrected count matrix for 14577 genes
## Calculating gene attributes
## Wall clock passed: Time difference of 17.54207 secs
## Determine variable features
## Centering data matrix
## Place corrected count matrix in counts slot
## Warning: Different cells and/or features from existing assay SCT
## Set default assay to SCT
## PC_ 1 
## Positive:  Umod, Egf, Tmem52b, Sult1d1, Sostdc1, Fabp3, Klk1, Prdx5, Cox6c, Foxq1 
##     Ly6a, Wfdc15b, Wfdc2, Slc25a5, Krt7, Cox5b, Ckb, Atp5g1, Cldn19, Cox4i1 
##     Ndufa4, Atp1a1, Ggt1, Cox8a, Atp5b, Chchd10, Cox6b1, Cox7a2, Uqcrq, Uqcr11 
## Negative:  Pappa2, Zfand5, Aard, Robo2, Wwc2, Pde10a, Mir6236, Itga4, Nadk2, Neat1 
##     Nos1, S100g, Col4a4, Ramp3, Sgms2, Col4a3, Ptgs2, Irx1, Tmem158, Itprid2 
##     Ranbp3l, Bmp3, Camk2d, Sdc4, Etnk1, Cdkn1c, Zbtb20, Srrm2, Hnrnpa2b1, Peg3 
## PC_ 2 
## Positive:  Fos, Junb, Jun, Hspa1b, Hspa1a, Btg2, Egr1, Zfp36, Atf3, Ier2 
##     Fosb, Socs3, Klf6, Klf2, Jund, Dnajb1, Dusp1, Tsc22d1, Rhob, Gadd45g 
##     H3f3b, Gadd45b, Ubc, Cebpd, Mt1, Mt2, Actb, H2bc4, H1f2, Klf4 
## Negative:  Mir6236, CT010467.1, Egf, Pappa2, Slc12a1, Umod, Etnk1, Wnk1, Nme7, Atp1b1 
##     Sfrp1, Robo2, Sptbn1, Col4a4, Sec14l1, Hsp90b1, Zbtb20, Dst, Itprid2, Tmem52b 
##     Lars2, Pde10a, Kcnq1ot1, App, Mal, Atrx, Col4a3, Tfcp2l1, Syne2, Nfe2l1 
## PC_ 3 
## Positive:  Egf, Mir6236, CT010467.1, Umod, Tmem52b, Neat1, Nme7, Fos, Jun, Malat1 
##     Slc12a1, Junb, Kcnq1ot1, Lars2, Egr1, Sult1d1, Foxq1, Etnk1, Dst, Wnk1 
##     Slc5a3, Atp1b1, Atf3, Fosb, Zfp36, Btg2, Atrx, Klf6, Pnisr, Ivns1abp 
## Negative:  Fth1, Ubb, Ftl1, Fxyd2, Ldhb, Car15, Prdx1, Cd63, Cd9, Eif1 
##     Mpc2, Mgst1, Mt1, Tmem213, Bsg, Clu, Mdh1, Aldoa, Itm2b, Ramp3 
##     Selenow, Tmem59, S100a1, Spp1, Tmem176b, Tspo, Tmbim6, Tle5, Atpif1, Wfdc2 
## PC_ 4 
## Positive:  Pappa2, Aard, Tmem52b, Umod, Egf, Sult1d1, Tmem158, Foxq1, Mcub, Ptgs2 
##     Dctd, Wwc2, Iyd, Car15, Ramp3, Ptprz1, Cd9, Hsp90b1, S100g, Wnk1 
##     Pth1r, Cdkn1c, Defb42, Clu, Ctsc, Tmsb4x, Tagln2, Bmp2, Col4a3, Tsc22d1 
## Negative:  Mt1, Apoe, Mt2, Aebp1, Sostdc1, Gpx6, Fxyd2, Ptger3, Neat1, Ckb 
##     Fgf9, Egfl6, Ivns1abp, Mfsd4a, Defb1, Car4, CT010467.1, Plet1, Fkbp11, Atp5md 
##     Chchd10, Cox6c, Atp5k, Igfbp5, Tmem213, Mgst3, Chka, Atp5mpl, Avpr1a, Lpl 
## PC_ 5 
## Positive:  CT010467.1, Hspa1b, Hspa1a, Klk1, Pappa2, Car15, Cldn10, Mir6236, Fth1, Lars2 
##     Jun, Hspa8, Wfdc2, Itm2b, Ftl1, Fau, Uba52-ps, Ly6a, Sfrp1, Aard 
##     Pik3r1, Mal, Eef1a1, Atp1a1, Ptger3, Atp1b1, Tspo, Hspb1, Tmem176a, Gpx4 
## Negative:  S100g, Actb, Tmem52b, Sdc4, Abhd2, Egr1, Uroc1, Tmsb10, Ppia, Atf3 
##     Serf2, Slc39a1, Alkbh5, Ndufb1-ps, Cebpd, Igfbp7, Cox6c, Atp5md, Ramp3, Ndufa3 
##     Wfdc15b, Kdm6b, Atp5k, Atp5e, Rbm47, Ldhb-ps, Ranbp3l, Ppm1h, Atp5mpl, Rcan1
## Computing nearest neighbor graph
## Computing SNN
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 11448
## Number of edges: 356414
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7087
## Number of communities: 26
## Elapsed time: 0 seconds
## 13:28:13 UMAP embedding parameters a = 0.9922 b = 1.112
## 13:28:13 Read 11448 rows and found 16 numeric columns
## 13:28:13 Using Annoy for neighbor search, n_neighbors = 30
## 13:28:13 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 13:28:13 Writing NN index file to temp file /var/folders/7n/x74qctp91rng390gx0z9hmd80000gn/T//RtmpIFjmQt/file54ff32e1dca1
## 13:28:13 Searching Annoy index using 1 thread, search_k = 3000
## 13:28:15 Annoy recall = 100%
## 13:28:15 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 13:28:16 Initializing from normalized Laplacian + noise (using RSpectra)
## 13:28:16 Commencing optimization for 200 epochs, with 483954 positive edges
## 13:28:16 Using rng type: pcg
## 13:28:18 Optimization finished
DimPlot(SO4)

DimPlot(SO4,split.by = "treatment")

SO4@meta.data$seurat_clusters <- factor(
  SO4@meta.data$seurat_clusters,
  levels = c("0", "1", "2", "3", "4", "6", "8", "9", "10", "12", "15", "22", "18", 
             "19", "21", "16", "13", "20", "23", "17", "14", "5", "7", "11", "24", "25")
)


Idents(SO4) <- SO4$seurat_clusters


## Simplifying into a few types

markers.to.plot2 <- c("Pappa2",
                      "Nos1",
                      "Egf",
                      "Umod",
                      "Fos",
                      "Junb",
                      "Cxcl10",
                      "Isg15",
                       "Sult1d1",
                      "Cldn19",
                      "Foxq1"
)

DimPlot(SO4)

DotPlot(SO4,
features = markers.to.plot2,
dot.scale = 8,
dot.min = 0,
scale = FALSE,
scale.max = 100,
scale.min = 0,
col.min = -2.5,
group.by = "seurat_clusters",
col.max = 2.5)+
coord_flip()

Idents(SO4) <- "seurat_clusters"

all_markers <- FindAllMarkers(SO4, only.pos = TRUE, min.pct = 0.5, logfc.threshold = 0.5)
## Calculating cluster 0
## Warning: `PackageCheck()` was deprecated in SeuratObject 5.0.0.
## ℹ Please use `rlang::check_installed()` instead.
## ℹ The deprecated feature was likely used in the Seurat package.
##   Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## For a (much!) faster implementation of the Wilcoxon Rank Sum Test,
## (default method for FindMarkers) please install the presto package
## --------------------------------------------
## install.packages('devtools')
## devtools::install_github('immunogenomics/presto')
## --------------------------------------------
## After installation of presto, Seurat will automatically use the more 
## efficient implementation (no further action necessary).
## This message will be shown once per session
## Calculating cluster 1
## Warning in FindMarkers.default(object = data.use, cells.1 = cells.1, cells.2 =
## cells.2, : No features pass logfc.threshold threshold; returning empty
## data.frame
## Calculating cluster 2
## Calculating cluster 3
## Calculating cluster 4
## Calculating cluster 6
## Calculating cluster 8
## Calculating cluster 9
## Calculating cluster 10
## Calculating cluster 12
## Calculating cluster 15
## Calculating cluster 22
## Calculating cluster 18
## Calculating cluster 19
## Calculating cluster 21
## Calculating cluster 16
## Calculating cluster 13
## Calculating cluster 20
## Calculating cluster 23
## Calculating cluster 17
## Calculating cluster 14
## Calculating cluster 5
## Calculating cluster 7
## Calculating cluster 11
## Calculating cluster 24
## Calculating cluster 25
all_markers %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1, p_val_adj < 0.05) %>%
    slice_head(n = 3) %>%
    ungroup() -> top3

DoHeatmap(SO4, features = top3$gene) + NoLegend()

5 Assigning Subclasses for clusters

SO4@meta.data <- SO4@meta.data %>%
  mutate(subclass_MD = dplyr::case_when(
    seurat_clusters %in% c(0,1,2,4,6,10,21,18,19,12,22) ~ "type_1a",
     seurat_clusters %in% c(3,9,15,21,8) ~ "type_1b",
    seurat_clusters %in% c(7,24,25,5,11) ~ "type_2a",
    seurat_clusters %in% c(14,17) ~ "type_2b",
      seurat_clusters %in% c(16,13) ~ "type_3a",
    seurat_clusters %in% c(20) ~ "type_3b",
    seurat_clusters == 23 ~ "type_4",
  ))


SO4@meta.data$subclass_MD <- factor(
  SO4@meta.data$subclass_MD, 
  levels = c("type_1a","type_1b", "type_2a", "type_2b","type_3a","type_3b","type_4")
)


SO4@meta.data <- SO4@meta.data %>%
  mutate(subclass2_MD = dplyr::case_when(
    seurat_clusters %in% c(0,1,2,3,4,6,8,9,10,21,15,22,18,19,21,12) ~ "type_1",
    seurat_clusters %in% c(7,24,25,5,11,14,17) ~ "type_2",
    seurat_clusters %in% c(16,13,20) ~ "type_3",
    seurat_clusters %in% c(23) ~ "type_4"
    
  ))

SO4@meta.data$subclass2_MD <- factor(
  SO4@meta.data$subclass2_MD, 
  levels = c("type_1", "type_2","type_3","type_4")
)

Idents(SO4) <- SO4@meta.data$subclass_MD

DimPlot(object = SO4, reduction = "umap", group.by = "subclass_MD", label = TRUE)

DimPlot(object = SO4, reduction = "umap", group.by = "subclass2_MD", label = TRUE)

DimPlot(SO4,group.by = "subclass2_MD",split.by = "treatment")

DimPlot(SO4,group.by = "seurat_clusters",split.by = "treatment")

library(plotly)

dim_plot <- DimPlot(SO4, reduction = "umap",group.by = "seurat_clusters") +
  ggtitle("UMAP Dimensional Reduction")
interactive_plot <- ggplotly(dim_plot)

# Display the interactive plot
interactive_plot
# Assign your markers as a vector
genes_to_plot <- c(
  "Fos", "Socs3", "Gadd45b", "Wee1", "Hspa1a", "Sat1", # TYPE 3
  "Cxcl12", "Itpr2", "Bmp4", "Casr", "Grin2c", "Irx3", "Rap1gap", "App", "Wwc1", # TYPE 2
  "Syt5", "Syn3", "Cacna1d", "Slc6a7", "Robo2", "Begain", # TYPE 1
  "Pappa1", "Nos1", "Ptgs2", "Bmp2", "Atp2a3" # SUBCLUSTERS within type 1
)

DoHeatmap(SO4, features = genes_to_plot, group.by = "seurat_clusters")
## Warning in DoHeatmap(SO4, features = genes_to_plot, group.by =
## "seurat_clusters"): The following features were omitted as they were not found
## in the scale.data slot for the SCT assay: Pappa1, Wwc1, Irx3, Grin2c

6 Viewing genes that differentiate each Type

markers.to.plot3 <- c("Pappa2",
                      "Egf",
                      "Umod",
                      "Fos",
                      "Junb",
                      "Cxcl10",
                      "Isg15"
)


DotPlot(SO4,
features = markers.to.plot3,
dot.scale = 8,
dot.min = 0,
scale = FALSE,
scale.max = 100,
scale.min = 0,
col.min = -2.5,
group.by = "subclass2_MD",
col.max = 2.5)+
coord_flip()

DotPlot(SO4,
features = markers.to.plot2,
dot.scale = 8,
dot.min = 0,
scale = FALSE,
scale.max = 100,
scale.min = 0,
col.min = -2.5,
group.by = "subclass_MD",
col.max = 2.5)+
coord_flip()

Idents(SO4) <- "seurat_clusters"

all_markers <- FindAllMarkers(SO4, only.pos = TRUE, min.pct = 0.5, logfc.threshold = 0.5)
## Calculating cluster 0
## Calculating cluster 1
## Warning in FindMarkers.default(object = data.use, cells.1 = cells.1, cells.2 =
## cells.2, : No features pass logfc.threshold threshold; returning empty
## data.frame
## Calculating cluster 2
## Calculating cluster 3
## Calculating cluster 4
## Calculating cluster 6
## Calculating cluster 8
## Calculating cluster 9
## Calculating cluster 10
## Calculating cluster 12
## Calculating cluster 15
## Calculating cluster 22
## Calculating cluster 18
## Calculating cluster 19
## Calculating cluster 21
## Calculating cluster 16
## Calculating cluster 13
## Calculating cluster 20
## Calculating cluster 23
## Calculating cluster 17
## Calculating cluster 14
## Calculating cluster 5
## Calculating cluster 7
## Calculating cluster 11
## Calculating cluster 24
## Calculating cluster 25
# Order by avg_log2FC

all_markers <- all_markers %>%

  arrange(desc(avg_log2FC)) %>%

  select(gene, everything())

 

# Split by cluster (ident column)

marker_list <- split(all_markers, all_markers$cluster)

 

# Sort names alphanumerically

sorted_cluster_names <- names(marker_list)[order(as.numeric(names(marker_list)))]

 

# Create a workbook

wb <- createWorkbook()

 

# Add each cluster as a new worksheet

for (cluster_name in sorted_cluster_names) {

  addWorksheet(wb, sheetName = cluster_name)

  writeData(wb, sheet = cluster_name, marker_list[[cluster_name]])

}

 

date <- format(Sys.Date(), "%Y%m%d")

 

# Save workbook

saveWorkbook(wb, here("jk_code", paste0(date, "_", "_SO4_FindAllMarkers_By_Cluster.xlsx")), overwrite = TRUE)

# Saving Seurat Object

save(SO4, file = here("jk_code", "JK_clean_MD.rds"))
all_markers %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1, p_val_adj < 0.05) %>%
    slice_head(n = 3) %>%
    ungroup() -> top3

DoHeatmap(SO4, features = top3$gene) + NoLegend()
## Warning in DoHeatmap(SO4, features = top3$gene): The following features were
## omitted as they were not found in the scale.data slot for the SCT assay: Nubp2